Structure and Diffusion of Nanoparticle Monolayers Floating at Liquid/ Vapor 

Interfaces: A Molecular Dynamics Study 



(N 

T 1 

O 

(N 



Shengfeng Chen^ and Gary S. Grest 
Sandia National Laboratories, Albuquerque, NM 87185, USA 
(Dated: May 2, 2012) 

Large-scale molecular dynamics simulations are used to simulate a layer of nanoparticles diffusing 
on the surface of a liquid. Both a low viscosity liquid, represented by Lennard- Jones monomers, 
and a high viscosity liquid, represented by linear homopolymers, are studied. The organization 
and diffusion of the nanoparticles are analyzed as the nanoparticle density and the contact angle 
between the nanoparticles and liquid are varied. When the interaction between the nanoparticles 
and liquid is reduced the contact angle increases and the nanoparticles ride higher on the liquid 
surface, which enables them to diffuse faster. In this case the short range order is also reduced 
as seen in the pair correlation function. For the polymeric liquids, the out-of-layer fluctuation is 
suppressed and the short range order is slightly enhanced. However, the diffusion becomes much 
slower and the mean square displacement even shows sub-linear time dependence at large times. The 
relation between diffusion coefficient and viscosity is found to deviate from that in bulk diffusion. 
Results are compared to simulations of the identical nanoparticles in 2-dimensions. 
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I. INTRODUCTION 

Nanoparticles at a liquid/vapor or liquid/liquid inter- 
face have attracted extensive attention during the past 
two decades.^ One motivation is that nanoparticles ad- 
sorbed at interfaces are found to be able to stabilize 
emulsions and foams<^"— Nanoparticles also self-assemble 
into various structures at an interface, which provides 
an efficient route to produce superlattices of nanopar- 
ticles of technological importancei^""— The advantage of 
this technique is that the assembly process can be fast 
and the resulting arrays are usually highly ordered. Due 
to their small size the adsorption energy of nanoparticles 
at an interface is typically only a few k^T , where /cb is 
the Boltzmann constant and T is the temperaturei^i^ 
The small adsorption energy implies that the adsorbed 
nanoparticles are highly dynamic and the self-assembly 
is quite reversible, the latter of which has been found to 
be crucial to form highly ordered arraysjS*^. Nanoparti- 
cles straddling an interface can also be regarded as living 
in a quasi 2-dimensional (2D) environment, which pro- 
vides a model system to study interesting problems such 
as phase transitions of 2D fluidsjiSrJ^ 

Many experimental studies have been devoted to inves- 
tigate the factors controlling the behavior of nanoparti- 
cles at an interface, including their size, surface morphol- 
ogy, shape, materials polarizability, and coatingsJ^"— 
Some of these factors influence the location and orien- 
tation of individual nanoparticles; others influence their 
mutual interactions and assembly geometry. However, 
in experiments these factors are usually intertwined to 
yield collective effects and it is difficult to single out the 
effect of each factor alone. This aspect is where molecu- 
lar dynamics (MD) simulations can play a useful role by 
studying the effect of one factor at a time to help elu- 
cidate experimental observations and uncover new phys- 
ical insights ji^Ti^ For example, Bresme et al. showed 
that Young's equation can be used to describe force bal- 
ance at nanoscale interfaces, but in certain cases line ten- 



sion should also be includedji^ Fenwick and Powell et al. 
showed that contrary to expectations, the collapse pres- 
sure measured in a typical Langmuir trough experiment 
should be independent of the contact anglei^i^ Recently, 
Cheung showed the importance of nanoparticle-liquid in- 
teractions and capillary waves in determining the stabil- 
ity of nanoparticles at liquid interfaces:^ 

Although MD simulations to date have revealed many 
important aspects of the physical behavior of nanoparti- 
cles at an interface, some aspects are still unclear. Par- 
ticularly, it is not clear how the structure and dynam- 
ics depend on contact angle and nanoparticle density in 
the low coverage regime. These behaviors are important 
since low density clusters can occur at the earlier stage 
of assembly and affect the morphology of the final dense 
layer 1^ The effect of liquid viscosity on interfacial diffu- 
sion is also not well understood and has only been studied 
recently.=2ii2£ In this paper we use MD to study the struc- 
ture and dynamics of a layer of nanoparticles floating at 
a liquid/vapor interface. We focus on the effect of vary- 
ing the contact angle, which is controlled by interactions 
between the nanoparticles and liquid, the nanoparticle 
density, and the liquid viscosity, respectively. Snapshots 
of some of the systems we have simulated are shown in 
Fig. [T] (only a part of the full simulation box is shown 
in each case). In Fig. [TJa) the contact angle between 
the nanoparticle and liquid is Oc = 137°. In this case 
the nanoparticles only slightly dip into the liquid which 
is composed of Lennard- Jones (LJ) monomers and is in 
equilibrium with its vapor phase. In Fig. [IJb) 9c = 29° 
and the nanoparticles are almost immersed in the liquid. 
In Fig. [IJc) 6c = 93° and the nanoparticles straddle the 
surface of the liquid, in this case composed of flexible 
linear polymer chains. Note that the LJ monomer liq- 
uid has a high vapor density, but the vapor density of 
the polymeric liquid is essentially 0.'^*' Also note that the 
polymeric liquid has a higher bulk density than that of 
the monatomic liquid and the liquid fllm in Fig. [IJc) is 
thinner than those in Fig.[TJa) and (b) since in our simu- 
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lations all liquids contain a similar number of monomers. 





FIG. 1: (Color online) Snapshots of nanoparticles floating at 
liquid/vapor interfaces: (a) monatomic LJ liquid and 9c = 
137°; (a) monatomic LJ liquid and 6^ = 29°; (c) 100-bead 
chain polymeric liquid and 6c = 93°. Only a small portion 
(lOOcr X 100(7 X lOOcr) of the simulation cell is shown in each 
snapshot. 



II. SIMULATION METHODOLOGY 

We placed a layer of nanoparticles at a liquid/ vapor 
interface as shown in Fig. [TJ Three liquid systems con- 
sisting of either LJ monomers or flexible linear chains of 
N LJ beads for = 10 and 100 were studied. In all 
three cases, the beads interact with each other through 
the standard LJ 12-6 potential 
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(1) 

where r is the distance between two beads, e is the unit of 
energy, and a is the diameter of beads. The interaction 
is truncated at = 3. Oct. For the two polymeric liquids, 
beads in the chain are connected by an additional finite 
extensible nonlinear elastic potential with a spring con- 
stant k — 30e/CT^ and maximum extent Rq — 1.5ct.^^ The 
liquid density p and shear viscosity r] for the L J monomer 
liquid IS p = QMa'^ and = l.Ol ± 0.03m/TCT,3a while 
for the polymer liquids, p = 0.89ct^"^ and r/ ~ 12±lm/rCT 
for N = 10, and p = QMa-^ and 77 = 180 ± IOto/tct for 
N = 100 ,33-34 

The nanoparticles are assumed to consist of a uniform 
distribution of atoms interacting with a L J potential. For 
spherical particles, their mutual interaction can then be 
determined analytically by integrating over all the inter- 
acting LJ atom pairs between the two particles i^^'^^ For 
nanoparticles with radii a, the interaction energy is given 
by 
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Here r is the center-to-center distance between two 
nanoparticles. The Hamaker constant A^^ — 
47r^ennP^o'^, where enn is the interaction strength be- 
tween the LJ atoms that make up the nanoparticles, and 



CT„ is the diameter and p„ the density of L J atoms in the 
nanoparticles. To reduce the number of parameters, we 
take e„n — e, ct„ = ct and p„ = I.Oct^"^, in which case 
^nn = 39.48e. In this paper we set a — IOct. 

The interaction between the LJ beads and nanoparti- 
cles is determined by integrating the interaction between 
a LJ bead and the LJ atoms within a nanoparticle, and 
the interaction potential Uns{r) is given by 
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where r is the center-to-center distance between the bead 
and nanoparticle, and the Hamaker constant A^s = 
247rensp„CT^j = 247rens for ct„ = ct and p„ = I.Oct^"^. 

Depending on the values of the Hamaker constant 
Ann and A^s, the nanoparticles can either be dispersed 
in the liquid or phase separate. As we are interested 
in studying nanoparticles that are hard-sphere like at 
the liquid/vapor surface, we truncate the nanoparticle- 
nanoparticle interaction so that it is purely repulsive. 
For a = IOct, this gives a cutoff = 20.427ct for J7nn('^)- 
Physically this corresponds to adding a short surfactant 
coating to the nanoparticles to avoid flocculation<^i^ 
For the interaction between the nanoparticles and LJ 
beads making up the liquid and vapor, we set Tc = 
a -t- 4ct = 14ct. The remaining free parameter A^s con- 
trols the solubility of the nanoparticles in the liquid. We 
choose such that nanoparticles phase separate to the 
liquid/vapor interface, in which case controls the con- 
tact angle 9c of the nanoparticles on the liquid surface. 
Results for 6c as a function of A^s are shown in Fig. [2] 
The data show that dc 180° as 0. However, 

when Ans exceeds certain critical value, which is approx- 
imately 85e for LJ monomers and 120e (130e) for 10-bead 
(100-bead) chains, 9c goes to and the nanoparticles dif- 
fuse into the liquid. In the intermediate range around 
90°, 9c decreases roughly linearly as increases. For a 
fixed ^ns, Sc increases as the chain length increases. This 
trend can be qualitatively understood as the result of less 
entropy gain when mixing the nanoparticles with longer 
chains. Correspondingly, a larger A^s, which represents 
the enthalpy contribution of solvation, is required to dis- 
perse the nanoparticles into a liquid of longer chains <^ 
The effects of and chain length on the behavior of the 
nanoparticles at the liquid/vapor interface are visualized 
in Fig. [11 

All MD simulations were performed using the 
LAMMPS simulation package The simulation cell 
is a rectangular box of dimensions x Ly x L^, where 
= 390.577CT, Ly = 451.0ct, and = 120ct for the 
LJ monomer liquid and 70ct for the polymeric liquids. 
The liquid/vapor interface is parallel to the x-y plane, 
in which periodic boundary conditions were employed. 
In the z direction, the LJ atoms and nanoparticles are 
confined by two flat walls at z = and z = Lz, re- 
spectively. Since the 2D packing of hard spheres is very 
sensitive to the aspect ratio of the enclosing box, we set 
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FIG. 2: (Color online) Contact angle 6c vs. nanoparticle- 
liquid interaction strength Ans for different liquids: 
monatomic LJ liquid (circles), 10-bead chain polymeric liquid 
(triangles), and 100-bead chain polymeric liquid (squares). 
Uncertainties in 6c are ~ 1° — 2°, smaller than the symbol 
size. Lines are guides to the eye. 



Lx/Ly = \/3/2 so that the hexagonal close-packing is 
favored. Each system contains more than 6 million LJ 
atoms to form a liquid layer with a thickness ^ SOcr 
and in equilibrium with its vapor phase, which in the LJ 
monomer case has a thickness ~ TOcr for Lz = 120a. For 
the polymeric liquids, the vapor density is and Lz was 
reduced accordingly. In order to investigate the effect of 
nanoparticle coverage, three systems were simulated with 
Np — 200, 240, and 320 nanoparticles, corresponding to 
2D density ^ = Np{2af / [L^Ly) = 0.45, 0.54, and 0.73, 
respectively. All these densities are well below the critical 
density 0.89 at the hquid/hexatic transition and 0.92 at 
the hexatic/solid transition of 2D hard disk fluids 42ri^ 

The LJ atoms interact with both upper and lower walls 
through a LJ 9-3 potential, which depends only on their 
distance z from the wall. 



for larger values of where the nanoparticles are dis- 
persed in the liquid.— 

The equations of motion were integrated using a 
velocity- Verlet algorithm with a time step 5t = 0.005r, 
where r = (T(m/e)^/^ and m is the mass of a LJ monomer. 
The nanoparticle with a radius a = IOct has a mass 
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4188. 79to. During the equilibration, the 



temperature T was held at LOe/fce by a Langevin ther- 
mostat weakly coupled to all LJ atoms with a damping 
constant T = 0.1r~^. Once the liquid/vapor interface 
was equilibrated, the Langevin thermostat was removed 
except for those liquid atoms within lOcr of the lower wall 
at = 0. Since the thickness of the liquid layer is ^ 50a 
and all nanoparticles are floating at the liquid/vapor in- 
terface in our simulations, their motion is not affected by 
the thermostat. 

For comparison, we also conducted MD simulations of 
5000 nanoparticles in a 2D box with L^jLy = 1. The size 
of the 2D box was varied to ensure the same density as for 
the nanoparticles at liquid/vapor interfaces. A Langevin 
thermostat with a damping constant F = O.lr"^ was 
used to keep the temperature at T = LOe/fce- The ther- 
mostat works as an implicit solvent. The interaction be- 
tween the nanoparticles is still given by Eq. ^ since 
capillary interactions are negligible for nanoparticles be- 
cause of the irrelevance of gravity.— We confirmed this 
treatment in our simulations by directly calculating and 
visualizing the liquid/vapor interface, which is flat (ex- 
cept for temporary capillary fluctuations) all the way to 
the contact line on the nanoparticle surface and shows 
no distortion at all with the presence of nanoparticles. It 
should be pointed out that for larger particles, capillary 
interactions can be taken into account even in 2D simula- 
tions by adding an effective capillary attraction between 
particles 4^ 

III. RESULTS AND DISCUSSION 
A. Out-of-Plane Fluctuations 



where ew = 
truncated at 
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2e. For the lower wall the interaction is 
Zc = 3.0(7, while at the upper wall is purely 
repulsive with z^ ~ 0.71476(7. Though all nanoparticles 
are confined to the liquid/ vapor interface and far from 
the two walls in most simulations, we also included a 
nanoparticle- wall potential of the form 



U{z) - A 



7560 
laz 



( 7a — z 

rln 



7a-\-z 
(z+aY 



z-\-a 



(5) 



where z is the distance of the center of a nanoparticle 
from the wall and = 144e. At both walls the po- 
tential is truncated at z = 10.57187(7 to make it purely 
repulsive for the nanoparticles. This potential is useful 



At a finite contact angle 9c, the nanoparticles strad- 
dle the liquid/ vapor interface and form a layer that is 
essentially 2D, with small thermal fluctuations in the 
z-direction normal to the interface. Direct calculation 
of the magnitude of fluctuations shows that it is almost 
independent of Oc and the nanoparticle density 4>, but 
strongly depends on the properties of the supporting liq- 
uid film, including its density, viscosity, and surface ten- 
sion. Because of thermal fiuctuations, the layer thickness 
is broadened and the individual nanoparticle height de- 
viates from the mean value of all nanoparticles. Such 
deviations are apparent in snapshots shown in Fig. [T] 
More quantitative results are shown in Fig. [31 where the 
probability density distribution P{5zp) is calculated as 
a function of the deviation 5zp from the instantaneous 
mean height of all nanoparticles. Note that the mean 
height itself fluctuates with time, but such fluctuations 
are excluded in Fig. [31 Including them would make the 
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FIG. 3: (Color online) Probability density distribution P{5zp) 
of nanoparticle positions in the z-direction, where 5zp is 
the deviation of position from the instantaneous mean of all 
nanoparticles. Data are for 4> = 0.54, dc ~ 93°, and various 
liquids: monatomic LJ liquid (circles), 10-bead chain poly- 
meric liquid (triangles), and 100-bead chain polymeric liquid 
(squares). Lines are the corresponding Gaussian fits. 

distributions even wider. The distributions in Fig. [3] all 
have a Gaussian shape, with a variance (1.75±0.08)cr for 
the monatomic LJ liquid, and ^ 0.7ct for the two poly- 
meric liquids. Since the nanoparticle diameter is 20ct, the 
height variations amongst the nanoparticles are less than 
10% of their diameter, which confirms that the nanopar- 
ticle layer is close to a 2D system. 

B. In-Plane Structure 

To characterize the structure in the plane of the 
nanoparticle layer, we computed the 2D radial distribu- 
tion function g{r) and the structure factor 5'(q), where 
q = 9a: By is a 2D wave- vector and tx and Cy are unit 

vectors along x and y directions, respectively. While g(r) 
and 5'(q) are related through a Fourier transform, for a fi- 
nite system it is easier to calculate each one directly. The 
calculation of g{r) involves counting the number of pairs 
of nanoparticles separated by distance r and is straight- 
forward. The structure factor is given by 

Siq) ^ Np^^eyi]i{ic^-v^n), (6) 

m.n 

where the sum is taken over all nanoparticle pairs indexed 
by m and n and separated by Tmn = — Tn • 

Results of g{r) are shown in Fig. 2] Data in Fig. HJa) 
are for different (/)'s at Oc — 93°. These results show the 
expected increase in the height of the first peak in g(r) 
as increases. The effect of 6c on g(r) is illustrated in 
Fig. IUb)-(d). All results consistently show that when 
9c is reduced, the locations of the peaks of g{r) move 
to larger r, which indicates that the nanoparticle layer 
is slightly denser at larger 9c- This trend is consistent 
with experimental results on the silica nanoparticle layer 
at the water-air interfaceJ^ It can be understood from 



a simple physical picture. At large 9c^ the nanoparticles 
ride high at the liquid/vapor interface and the separation 
between nanoparticles is solely controlled by the hard 
core repulsion. However, at small 0c, the nanoparticles 
are partially coated by the liquid, which increases their 
effective size. So as 9c decreases the mean separation 
between two nanoparticles increases. 
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FIG. 4: (color online) The pair distribution function g{r). 
Lines in all main panels are for the monatomic LJ liquid: 
(a) Be = 93° and (j) = 0.45 (solid), 0.54 (dashed), 0.73 (dot- 
ted); (b) (j) = 0.73 and Oc = 137° (solid), 93° (dashed), 51° 
(dotted); (c) <j) = 0.54 and 6^ = 137° (solid), 93° (dashed), 
51° (dotted); (d) </> = 0.45 and 6^ = 137° (solid), 93° 
(dashed), 51° (dotted). Inset of (b): The dashed line is for 
the monatomic LJ liquid and the dash-dotted line is for the 
polymeric liquid consisting of 100-bead chains; for both lines 
(j> = 0.54 and 6^ = 93°; 

Figure [DJb) shows that at (/) — 0.73 the peaks of g{r) 
grow as 9c decreases, indicating that the nanoparticle 
layer becomes more ordered at smaller 9c. The data in 
Fig. Etc) and (d) for 9c = 137° and 93° also show this 
trend. However, the first peak of g{r) becomes lower 
when 9c is further reduced to 51°, though other peaks 
become higher. The apparent reduction in the first peak 
of g(r) at 0c = 51° is due to a finite size effect. The cal- 
culation of S'(q) indicates that all peaks grow as 9c de- 
creases and confirms that the nanoparticle layer exhibits 
stronger local order at smaller 9c- 

As 9c continues to decrease to 0, we have observed the 
loss of nanoparticles at the interface as some of them dif- 
fuse into the liquid, though 9c for an individual nanopar- 
ticle may still be finite. This effect is more significant 
at higher nanoparticle density. For example, at 9c = 29° 
more than 18% nanoparticles were dispersed in the liquid 
at the end of MD runs for cj) = 0.73; while for (p = 0.45, 
about 8% have diffused into the liquid. The loss of 
nanoparticles at the interface for a finite 9c is a conse- 
quence of their small activation energy, which is typically 
at the order of several k-QTX 
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Result of g{r) for the polymeric liquid consisting of 
100-bead chains is included in the inset of Fig. IH^c), to- 
gether with the result for the monatomic LJ liquid at 
the same dc and (p. The peaks of g{r) are slightly higher 
for the polymeric liquid than those for the monatomic 
LJ liquid, indicating that local order is slightly stronger 
in the former case. As shown in Fig. [31 the out-of-layer 
fluctuation is much smaller for the polymeric liquids, i.e., 
the nanoparticle layer is more 2D-like. As a consequence, 
the nanoparticle density is effectively higher, which leads 
to higher peaks in g{r) for the polymeric liquids. 

A comparison of g{r) for nanoparticles at the liq- 
uid/vapor interface and in 2D with an implicit solvent 
is shown in Fig. [5] for three densities. The peaks in g{r) 
for nanoparticles at the interface are clearly higher and 
decay slower than those for nanoparticles in 2D, indicat- 
ing stronger local order in the former case. The locations 
of the peaks move towards larger r for nanoparticles at 
the interface. The shift is between 2 to 2.5a, and re- 
flects the liquid coating that is about 1 cr in thickness on 
each nanoparticle. The coating makes the nanoparticles 
effectively larger than the bare ones in 2D simulations. 
Therefore, the actual nanoparticle density is effectively 
higher at the interface, which leads to higher peaks in 
g{r). However, if the radius of nanoparticles in 2D was 
increased to reflect this coating, then we would expect 
higher peaks in g{r) at a given (/) from 2D simulations. 
This trend can be derived indirectly from the previous 
comparison of ^(r) between the monatomic LJ liquid and 
polymeric liquids. Higher peaks are found for the latter 
since the nanoparticle layer there is more 2D-like. 
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FIG. 5: (color online) The pair distribution function g{r) 
for nanoparticles at the liquid/ vapor interface for the LJ 
monatomic liquid at 6c = 93° (dashed lines) compared to 
results of 2D simulations (solid lines) for (j) — 0.45 (bottom), 
0.54 (middle), and 0.73 (top). 

Additional information on the in-plane structure of the 
nanoparticle layer can be obtained from 5'(q). A den- 
sity plot of 5'(q) in the Qx-Qy plane is shown in Fig. [6] 
for 9c = 93° and the monatomic LJ liquid at ^ = 0.54 
and 0.73, and the 100-bead chain polymeric liquid at 



4> = 0.54. At low density (j) = 0.45 and 0.54, the local 
structure is almost isotropic, indicating a fluid-like state 
of the floating layer. However, at = 0.73, which is still 
lower than the critical density for the fluid/solid tran- 
sition of 2D hard sphere systems, the hexagonal close- 
packing feature of the local structure is rather apparent 
as shown in the middle panel of Fig. |6l Smaller out- 
of-plane fluctuations of nanoparticles on the surface of 
the polymeric liquid lead to stronger local order in the 
nanoparticle layer at a given </> and 9c, indicated by higher 
peaks in 5(q). 
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FIG. 6: (color online) Density plots of ^(q) in the Qx-Qy plane 
for 8c = 93° and the monatomic LJ liquid at (j) — 0.54 (left) 
and 0.73 (middle), and the 100-bead chain polymeric liquid 
at (/> = 0.54 (right). 

C. Orientational Order 

A simple measure of orientational order of the float- 
ing nanoparticle layer is provided by the fraction of 
nanoparticles, /g, which have exactly six neighbors as 
in a close-packed hexagonal lattice. Here a simple cri- 
terion is adopted to identify nearest neighbors as those 
within a cut-off radius Vc of a given nanoparticle. We 
chose Tc = 35tT, roughly corresponding to the location of 
the first valley of g{r) as shown in Fig. 01 After nearest 
neighbors were found, we computed the Nelson-Halperin 
(N-H) order parameter, which can be expressed as 




where the first sum is over all Np nanoparticles, Uj is the 
number of nearest neighbors of the j-th nanoparticle, the 
second sum is over rij nearest neighbors, and 9jk is the 
angle formed by the bond between a nearest neighbor 
pair (j, k) and a fixed axis. For an ideal hexagonal lat- 
tice me — 1, and it decreases to zero as the local bond 
disorder increases. 

Results of /e and mg vs. (f) are shown in Fig. [71 for 
nanoparticles fioating on the surface of the monatomic 
LJ liquid at various contact angles. Results of 2D sim- 
ulations are also included. As expected, in all cases /e 
and mg increase with (j), indicating the development of 
orientational order. Figure [3 also shows that at a given 
(j), the layer is more locally ordered at smaller 9c as in- 
dicated by larger values of mg and /g. The increase in 
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FIG. 7: (a) The fraction of nanoparticles, /e, having exactly 
six neighbors and (b) the N-H order parameter mg vs. the 
nanoparticle density (f) for nanoparticles on the surface of the 
monatomic LJ liquid at various contact angles: 9c ~ 137° 
(circles); 9c = 93° (triangles); 9c — 51° (squares). Results of 
2D simulations are shown with diamonds. Lines are guides to 
the eye. 



for finding nearest neighbors needed for the calculation 
of the N-H order parameter, even when the nanoparticle 
density is low. 

D. Diffusion Coefficient 




t/T 



FIG. 8: (a) Mean square displacement (Ar^) vs. time t. The 
top three data sets are for the monatomic LJ liquid at (p = 
0.45 and 9c = 137°; rj) = 0.54 and 9c = 137°; (j) = 0.54 and 
9c ~ 93°. The second to bottom (bottommost) data set is for 
the polymeric liquid consisting of 10-bead (100-bead) chains 
at (A = 0.54 and 61^ = 93°. 



orientational order as 6c is reduced is more dramatic at 
higher (j). This trend persists in our simulations when 9c 
is reduced even further (e.g., 9c = 29°). In this case since 
some nanoparticles are eventually absorbed into the liq- 
uid at long times, we can only calculate fe and ttiq in the 
early stages of simulations when all nanoparticles are still 
at the liquid/vapor interface. In this early time regime 
systems follow the same trend that at a given (j), values of 
fe and me are larger and the order is stronger at smaller 
9c- The trend is consistent with that found from g{r) as 
shown in Fig. S) Data in Fig. [7] also show that orienta- 
tional order is stronger for nanoparticles at the interface 
than for those in 2D since the liquid coating in the former 
case makes nanoparticles effectively larger and their den- 
sity higher. The same reason also leads to enhancement 
of translational order as shown in Fig. [5l 

The N-H order parameter ttiq shown in Fig. [7] was cal- 
culated with nearest neighbors identified with a simple 
cut-off criterion. To verify that this criterion leads to 
an accurate estimate of orientational order, we also con- 
ducted a Voronoi analysis of the packing geometry of 
the nanoparticle layeri^ Nearest neighbors were identi- 
fied as those sharing common sides in the Voronoi con- 
struction. The N-H order parameter was then computed 
with Eq. ([7]). Results from the Voronoi analysis typically 
agree within a few percent with those from the simpler 
cut-off approach. Thus the cut-off criterion is adequate 



The in-plane motion of nanoparticles at the liq- 
uid/vapor interface is characterized by the mean square 
displacement (Ar^) = (Ax^ 4- Ay^), which is shown in 
Fig. [5] for various cases. At very short times, the motion 
of nanoparticles is nearly ballistic with (Ar^) = w^t^, 
where w is a velocity and t is time. At larger times, the 
motion becomes diffusive and (Ar^) = ADt, where D is 
a diffusion coefficient and 4 is the prefactor for 2D dif- 
fusion. Most of our results fit to this classical picture, 
as shown in Fig. [S] In all cases, v ~ 0.02a /t and is 
nearly independent of 9c and (f). The intersection of the 
ballistic and diffusive regime defines a ballistic time scale 
tb = AD/v^. Results for D and <{, are shown in Fig. [9] as 
a function of the true nanoparticle density cf) for 4 values 
of 9c- At small 9c (e.g., 9c = 29°), because of the loss of 
nanoparticle from the interface, (j) is lower than the nom- 
inal density (j) calculated from Np. In other cases (f> = (f)- 
Figure [5] shows that in all cases, D and tf, decrease ap- 
proximately linearly with the nanoparticle density and 
the slope is steeper for larger 9c- At the same 0, both 
D and ti, decrease as 9c decreases. This is due to the 
fact that the attraction between the nanoparticles and 
liquid is enhanced to make 9c smaller. As a consequence, 
nanoparticles become more immersed into the liquid and 
need to plow through more liquid in order to move, which 
makes their diffusion more difficult. At = 0.45, 13 and 
tb decrease by a factor of 3 when 9c is reduced from 137° 
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to 29°. At (j) = 0.73, the reduction is almost 6-fold for 
the same change in 9c- 
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FIG. 9: (a) Diffusion coefficient D and (b) ballistic time ti, 
vs. the true nanoparticle density (p at various contact angles: 
9c = 137° (circles), 6c = 93° (triangles), Oc = 51° (squares), 
9c = 29° (diamonds). <j> is the same as except for 9c = 29°, 
where (j) < (j). Lines are guides to the eye. 

For the 10-bead chain polymeric liquid, (Ar^) grows 
linearly with at early times and shows diffusive behav- 
ior at large times. However the diffusion coefficient is 
reduced significantly compared to the LJ monomer case. 
For example, at (/> = 0.54 and 9c = 93°, the value of D 
is ~ 1 X lO^^cr^/r for the 10-bead polymeric liquid com- 
pared to ^ 7 X 10~^(t^/t for the monatomic LJ liquid. 
Note that Figl2] shows that the interaction strength be- 
tween the nanoparticles and liquid, dictated by A„s, has 
to be adjusted to ensure the same 9c for the polymeric 
and monatomic LJ liquid. For example, 9c = 93° for the 
monatomic LJ liquid at A^s = 60e/a^, while we need to 
increase A^s to 79e/cr^ to get the same 9c for the 10-bead 
polymeric liquid. If A^s is held fixed at A„s = 60e/(T^ for 
the 10-bead polymeric liquid, then 9c increases to 113°. 
In this case _D is ^ 1.6 x \Q^^a^ /t. These data are 
summarized in Table ID If we compare the values of D 
between the monatomic LJ and 10-bead polymeric liq- 
uid at A-ns — 60e/(T^, i.e., at the same nanoparticle-bead 
interaction, we find that D is reduced by a factor about 
4. Note that the viscosity r] of the 10-bead polymeric 
liquid is approximately 12 times larger than that of the 
monatomic LJ liquid. If we assume a scaling relation 
between D and rj as D ^ then a ~ 0.6 is much 

less than 1, which is in contrast with the situation of 
free bulk Brownian diffusion where a = 1 is expected. 



liquid Ans (e/g-') Oc (°) n [m/aT) D (a'/r) 

monatomic 60 93 1.01 7.0 x 10~^ 

10-bead 60 113 12 1.6 x 10"^ 

10-bead 79 93 12 1.0 x 10"^ 



TABLE I: Comparison of values of contact angle (9c), vis- 
cosity [rf], and diffusion coefficient [D) at = 0.54 for the 
monatomic LJ and 10-bead polymeric liquid. 

However, our results are consistent with a recent experi- 
mental measurement of the diffusion coefficient of single 
nanoparticles at water-oil interfaces, where a — 0.44 was 
found 1^ 

For the 100-bead chains, (Ar^) ^ P at early times, and 
then crosses over to a sub-diffusive regime as shown in 
Fig.m This deviation from the expected long-time linear 
dependence on time continued out to the longest times 
we are presently able to simulate. In this sub-diffusive 
regime, the mean square displacement for the 100-bead 
chains can be better fit to (Ar^) ~ as shown in 

Fig. [51 The reason underlying this sub-diffusive regime 
is unclear. 

V. CONCLUSIONS 

In this paper we studied nanoparticles fioating at liq- 
uid/vapor interfaces with MD simulations. Both the 
low-viscous monatomic LJ liquid and the high-viscous 
polymeric liquids composed of flexible linear chains were 
studied. We showed that as the attraction between the 
nanoparticles and liquid is increased, the contact angle 
is reduced and the nanoparticles are more wetted by the 
liquid. At the same time the short range order of the 
nanoparticle layer is slightly enhanced and the interfa- 
cial diffusion coefficient of the nanoparticles is greatly 
reduced. Our results further showed that both the trans- 
lational and orientational order of the nanoparticle layer 
grow quickly and the nanoparticle diffusion slows down 
dramatically as the nanoparticle density is increased. 
Comparisons with results of 2D simulations revealed that 
the main effect of liquid on the nanoparticles is to provide 
a coating which makes their effective size larger than that 
of the bare ones. Otherwise, the nanoparticle layer at the 
liquid/vapor interface is close to a 2D system, though the 
out-of-planc fluctuations can be as large as 10% of their 
diameter. The simulations with more viscous polymeric 
liquids showed that the the out-of-plane fluctuations of 
the nanoparticles are strongly suppressed even for rel- 
atively short chains (e.g., 10-bead chains). The local, 
short range order is slightly enhanced at a given con- 
tact angle and nanoparticle density as the chain length 
increases. At the same time, the nanoparticle diffusion 
becomes slower for more viscous liquids and even shows 
sub-diffusive behavior at large times for highly viscous 
liquids. The diffusion coefficient scales inversely with the 
viscosity however with an exponent less than 1. 
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